Capitulo 4

10

Debemos de utilizar la base de datos Weekly, esta base de datos contiene información del rendimiento porcentual semanal para el índice bursátil S&P 500 entre 1990 y 2010.

# Se carga la base de datoa a utilizar
library(ISLR)
# Estadístico básicos
summary(Weekly)
##       Year           Lag1               Lag2               Lag3         
##  Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
##  1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580  
##  Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410  
##  Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472  
##  3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
##  Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
##       Lag4               Lag5              Volume       
##  Min.   :-18.1950   Min.   :-18.1950   Min.   :0.08747  
##  1st Qu.: -1.1580   1st Qu.: -1.1660   1st Qu.:0.33202  
##  Median :  0.2380   Median :  0.2340   Median :1.00268  
##  Mean   :  0.1458   Mean   :  0.1399   Mean   :1.57462  
##  3rd Qu.:  1.4090   3rd Qu.:  1.4050   3rd Qu.:2.05373  
##  Max.   : 12.0260   Max.   : 12.0260   Max.   :9.32821  
##      Today          Direction 
##  Min.   :-18.1950   Down:484  
##  1st Qu.: -1.1540   Up  :605  
##  Median :  0.2410             
##  Mean   :  0.1499             
##  3rd Qu.:  1.4050             
##  Max.   : 12.0260

Tenemos una base de datos con 1089 observaciones sobre las siguientes 9 variables.

  • Year: El año en que se registró la observación.
  • Lag1: Porcentaje de retorno de la semana anterior.
  • Lag2: Porcentaje de retorno de 2 semanas anteriores
  • Lag3: de retorno de 3 semanas anteriores
  • Lag4: Porcentaje de retorno de 4 semanas anteriores
  • Lag5: Porcentaje de retorno de 5 semanas anteriores
  • Volume: Volumen de acciones negociadas (número promedio de acciones diarias negociadas en miles de millones)
  • Today: Porcentaje de retorno para esta semana
  • Direction: Un factor con niveles hacia abajo y hacia arriba que indica si el mercado tuvo un rendimiento positivo o negativo en una semana determinada
library(psych)
# Gráfico de dispersión
pairs.panels(Weekly[,-9], 
             method = "pearson",
             hist.col = "#00AFBB",
             density = TRUE,  
             )

Como podemos apreciar, en general no pareciera que hubiesen muchas variables correlacionadas, pero si podemos apreciar una correlación positiva entre el año y la variable Volume. La variable volume hace referencia al volumen de acciones negociadas (número promedio de acciones diarias negociadas en miles de millones), por lo que a medida que pasa el tiempo, podemos apreciar como se aumenta también el volumen de acciones negociadas.

Le ejecutamos a todos los datos, utilizando Direction como la variable respuesta y las variables lag más la variable volume como predictoras, una regresión logística

weeklyrlo <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Weekly, family = "binomial") 
summary(weeklyrlo)
## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = "binomial", data = Weekly)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6949  -1.2565   0.9913   1.0849   1.4579  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.26686    0.08593   3.106   0.0019 **
## Lag1        -0.04127    0.02641  -1.563   0.1181   
## Lag2         0.05844    0.02686   2.175   0.0296 * 
## Lag3        -0.01606    0.02666  -0.602   0.5469   
## Lag4        -0.02779    0.02646  -1.050   0.2937   
## Lag5        -0.01447    0.02638  -0.549   0.5833   
## Volume      -0.02274    0.03690  -0.616   0.5377   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1496.2  on 1088  degrees of freedom
## Residual deviance: 1486.4  on 1082  degrees of freedom
## AIC: 1500.4
## 
## Number of Fisher Scoring iterations: 4

Podemos apreciar que la variable lag2, es la única variables significativa.

library("caret")
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
predichos=as.factor(ifelse(test = weeklyrlo$fitted.values > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Weekly$Direction)
matriz<-confusionMatrix(predichos, verdaderos)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down  Up
##       Down   54  48
##       Up    430 557
##                                          
##                Accuracy : 0.5611         
##                  95% CI : (0.531, 0.5908)
##     No Information Rate : 0.5556         
##     P-Value [Acc > NIR] : 0.369          
##                                          
##                   Kappa : 0.035          
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.11157        
##             Specificity : 0.92066        
##          Pos Pred Value : 0.52941        
##          Neg Pred Value : 0.56434        
##              Prevalence : 0.44444        
##          Detection Rate : 0.04959        
##    Detection Prevalence : 0.09366        
##       Balanced Accuracy : 0.51612        
##                                          
##        'Positive' Class : Down           
## 
Reference
Predicted Down Up
Down A B
Up C D
  • Lo primero que podemos apreciar es que la precisión del modelo es de un 56.11% (Accuracy).
  • También podemos ver que la sensitividad (\(Sensitivity = \frac{A}{A + C}\)), es de solamente un 11.15%, lo cual quiere decir que en semanas en donde el mercado es bajista.
  • La especificidad (\(Specificity = \frac{D}{B + D}\)), es de un 92.06%, es la precisión a la hora de predecir semanas alcistas
attach(Weekly)
w_entrenamiento <- (Year <= 2008) 
weeklyrlo2 <- glm(Direction ~ Lag2, data = Weekly, family = "binomial", subset = w_entrenamiento)
summary(weeklyrlo2)
## 
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = Weekly, 
##     subset = w_entrenamiento)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.536  -1.264   1.021   1.091   1.368  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.20326    0.06428   3.162  0.00157 **
## Lag2         0.05810    0.02870   2.024  0.04298 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1354.7  on 984  degrees of freedom
## Residual deviance: 1350.5  on 983  degrees of freedom
## AIC: 1354.5
## 
## Number of Fisher Scoring iterations: 4
Weekly.20092010 <- Weekly[!w_entrenamiento, ]
Direction.20092010 <- Direction[!w_entrenamiento]

p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down Up
##       Down    9  5
##       Up     34 56
##                                          
##                Accuracy : 0.625          
##                  95% CI : (0.5247, 0.718)
##     No Information Rate : 0.5865         
##     P-Value [Acc > NIR] : 0.2439         
##                                          
##                   Kappa : 0.1414         
##                                          
##  Mcnemar's Test P-Value : 7.34e-06       
##                                          
##             Sensitivity : 0.20930        
##             Specificity : 0.91803        
##          Pos Pred Value : 0.64286        
##          Neg Pred Value : 0.62222        
##              Prevalence : 0.41346        
##          Detection Rate : 0.08654        
##    Detection Prevalence : 0.13462        
##       Balanced Accuracy : 0.56367        
##                                          
##        'Positive' Class : Down           
## 
  • Podemos apreciar que la precisión del modelo es de un 62.50%
  • Cuando la semana del mercado es bajista, la precisión es de un 20.93%
  • Cuando la semana del mercado es alcista, la precisión es de un 91.80%
library(MASS)
ldapre <- lda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
ldapre
## Call:
## lda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581
## 
## Coefficients of linear discriminants:
##            LD1
## Lag2 0.4414162
pred.lda <- predict(ldapre, Weekly.20092010)
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down Up
##       Down    9  5
##       Up     34 56
##                                          
##                Accuracy : 0.625          
##                  95% CI : (0.5247, 0.718)
##     No Information Rate : 0.5865         
##     P-Value [Acc > NIR] : 0.2439         
##                                          
##                   Kappa : 0.1414         
##                                          
##  Mcnemar's Test P-Value : 7.34e-06       
##                                          
##             Sensitivity : 0.20930        
##             Specificity : 0.91803        
##          Pos Pred Value : 0.64286        
##          Neg Pred Value : 0.62222        
##              Prevalence : 0.41346        
##          Detection Rate : 0.08654        
##    Detection Prevalence : 0.13462        
##       Balanced Accuracy : 0.56367        
##                                          
##        'Positive' Class : Down           
## 
  • La precisión de este modelo es de 62.50%
  • Cuando en la semana el mercado es bajista la precisión es de 20.93%
  • Cuando en la semana el mercado es alcista la precisión es de 91.80%
qdapre <- qda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
qdapre
## Call:
## qda(Direction ~ Lag2, data = Weekly, subset = w_entrenamiento)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down Up
##       Down    0  0
##       Up     43 61
##                                           
##                Accuracy : 0.5865          
##                  95% CI : (0.4858, 0.6823)
##     No Information Rate : 0.5865          
##     P-Value [Acc > NIR] : 0.5419          
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : 1.504e-10       
##                                           
##             Sensitivity : 0.0000          
##             Specificity : 1.0000          
##          Pos Pred Value :    NaN          
##          Neg Pred Value : 0.5865          
##              Prevalence : 0.4135          
##          Detection Rate : 0.0000          
##    Detection Prevalence : 0.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : Down            
## 
  • En este caso la predicción del modelo es de un 58.65%
  • La precisión cuando es mercado es bajista es de un 0%
  • La precisión cuando el mercado es alcista es de un 100%
library(class)
entrenamiento.x <- as.matrix(Lag2[w_entrenamiento])
prueba.x <- as.matrix(Lag2[!w_entrenamiento])
entrenamiento.Direction <- Direction[w_entrenamiento]
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 1)
confusionMatrix(pred.knn, Direction.20092010)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Down Up
##       Down   21 30
##       Up     22 31
##                                           
##                Accuracy : 0.5             
##                  95% CI : (0.4003, 0.5997)
##     No Information Rate : 0.5865          
##     P-Value [Acc > NIR] : 0.9700          
##                                           
##                   Kappa : -0.0033         
##                                           
##  Mcnemar's Test P-Value : 0.3317          
##                                           
##             Sensitivity : 0.4884          
##             Specificity : 0.5082          
##          Pos Pred Value : 0.4118          
##          Neg Pred Value : 0.5849          
##              Prevalence : 0.4135          
##          Detection Rate : 0.2019          
##    Detection Prevalence : 0.4904          
##       Balanced Accuracy : 0.4983          
##                                           
##        'Positive' Class : Down            
## 
  • La precisión del modelo es de un 50%
  • Cuando el mercado es bajista la precisión es de un 48.84%
  • Cuando el mercado es alcista la precisión es de un 50.82%
  1. Si analizamos las precisiones de los modelos obtuvimos los siguientes resultados:
GLN LDA QDA KNN
Precisión 62.5 62.5 58.65 50
Sensitividad 20.93 20.93 0 48.84
Especificidad 91.80 91.80 100 50.82

Por lo que podemos concluir que tanto el modelo GLN como LDA fueron los que mejor resultados dieron.

A continuación evaluamos KNN con los valores de K = 10, 20, 50 y 100

# Para K = 10
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 10)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table
##           Reference
## Prediction Down Up
##       Down   16 20
##       Up     27 41
## [1] "Accuracy: 0.5481 Sensitivity: 0.3721 Specificity: 0.6721"
# Para k = 20
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 20)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table
##           Reference
## Prediction Down Up
##       Down   23 21
##       Up     20 40
## [1] "Accuracy: 0.6058 Sensitivity: 0.5349 Specificity: 0.6557"
# Para k = 50
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 50)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table
##           Reference
## Prediction Down Up
##       Down   20 22
##       Up     23 39
## [1] "Accuracy: 0.5673 Sensitivity: 0.4651 Specificity: 0.6393"
# Para k = 50
pred.knn <- knn(entrenamiento.x, prueba.x, entrenamiento.Direction, k = 100)
z <-confusionMatrix(pred.knn, Direction.20092010)
z$table
##           Reference
## Prediction Down Up
##       Down   10 11
##       Up     33 50
## [1] "Accuracy: 0.5769 Sensitivity: 0.2326 Specificity: 0.8197"

Ahora probaremos una regresión logística con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume

# para volume
weeklyrlo2 <- glm(Direction ~ Volume, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   31 47
##       Up     12 14
## [1] "Accuracy: 0.4327 Sensitivity: 0.7209 Specificity: 0.2295"
# para lag2 + lag3
weeklyrlo2 <- glm(Direction ~ Lag2 + Lag3, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    8  4
##       Up     35 57
## [1] "Accuracy: 0.625 Sensitivity: 0.186 Specificity: 0.9344"
# para lag1 + lag2
weeklyrlo2 <- glm(Direction ~ Lag1 + Lag2, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    7  8
##       Up     36 53
## [1] "Accuracy: 0.5769 Sensitivity: 0.1628 Specificity: 0.8689"
# para lag1 + lag2 + volume
weeklyrlo2 <- glm(Direction ~ Lag1 +  Lag2 + Volume, data = Weekly, family = "binomial", subset = w_entrenamiento)
p <- predict(weeklyrlo2, Weekly.20092010, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = "Up", no = "Down"))
verdaderos=as.factor(Direction.20092010)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   27 33
##       Up     16 28
## [1] "Accuracy: 0.5288 Sensitivity: 0.6279 Specificity: 0.459"

Ahora probaremos una LDA con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume

# Para volume
ldapre <- lda(Direction ~ Volume, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   33 47
##       Up     10 14
## [1] "Accuracy: 0.4519 Sensitivity: 0.7674 Specificity: 0.2295"
# Para lag2 + lag3
ldapre <- lda(Direction ~ Lag2 + Lag3, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    8  4
##       Up     35 57
## [1] "Accuracy: 0.625 Sensitivity: 0.186 Specificity: 0.9344"
# Para lag1 + lag2
ldapre <- lda(Direction ~ Lag1 + Lag2, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    7  8
##       Up     36 53
## [1] "Accuracy: 0.5769 Sensitivity: 0.1628 Specificity: 0.8689"
# Para lag1 + lag2 + volume
ldapre <- lda(Direction ~ Lag1 + Lag2 + Volume, data = Weekly, subset = w_entrenamiento)
pred.lda <- predict(ldapre, Weekly.20092010)
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   27 33
##       Up     16 28
## [1] "Accuracy: 0.5288 Sensitivity: 0.6279 Specificity: 0.459"

Ahora probaremos una QDA con las diferentes configuraciones de variables predictoras: volume, lag2 + lag3, lag1 + lg2, lag1 + lag2 + volume

# para volume
qdapre <- qda(Direction ~ Volume, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   43 59
##       Up      0  2
## [1] "Accuracy: 0.4327 Sensitivity: 1 Specificity: 0.0328"
# para lag2 + lag3
qdapre <- qda(Direction ~ Lag2 + Lag3, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    4  2
##       Up     39 59
## [1] "Accuracy: 0.6058 Sensitivity: 0.093 Specificity: 0.9672"
# para lag1 + lag2
qdapre <- qda(Direction ~ Lag1 + Lag2, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down    7 10
##       Up     36 51
## [1] "Accuracy: 0.5577 Sensitivity: 0.1628 Specificity: 0.8361"
# para lag1 + lag2 + volume
qdapre <- qda(Direction ~ Lag1 + Lag2 + Volume, data = Weekly, subset = w_entrenamiento)
pred.qda <- predict(qdapre, Weekly.20092010)
matriz<-confusionMatrix(pred.qda$class, verdaderos)
matriz$table
##           Reference
## Prediction Down Up
##       Down   31 44
##       Up     12 17
## [1] "Accuracy: 0.4615 Sensitivity: 0.7209 Specificity: 0.2787"

11

Auto <- Auto[c(-10,-11,-12,-13,-14,-15,-16,-17,-18)]
mpg01 <- rep(0, length(Auto$mpg))
mpg01[Auto$mpg > median(Auto$mpg)] <- 1
Auto <- data.frame(Auto, mpg01)
str(Auto)
## 'data.frame':    392 obs. of  10 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
##  $ weight      : num  3504 3693 3436 3433 3449 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
##  $ mpg01       : num  0 0 0 0 0 0 0 0 0 0 ...
# Sabiendo que name es un factor, lo excluimos del análisis

pairs.panels(Auto[,-9], 
             method = "pearson",
             hist.col = "#00AFBB",
             density = TRUE,  
             )

Podemos apreciar que existe una correlación negativa con las variables cylinders, displacement, horsepower y weight

x <- sample(1:dim(Auto)[1], size=dim(Auto)[1]*0.75)
train <- Auto[x,]
test = Auto[-x,]
autos.lda = lda(mpg01~displacement+weight+cylinders+horsepower, data=train)
autos.lda
## Call:
## lda(mpg01 ~ displacement + weight + cylinders + horsepower, data = train)
## 
## Prior probabilities of groups:
##   0   1 
## 0.5 0.5 
## 
## Group means:
##   displacement   weight cylinders horsepower
## 0     271.4830 3601.537  6.789116  129.09524
## 1     118.0782 2352.027  4.217687   79.47619
## 
## Coefficients of linear discriminants:
##                        LD1
## displacement -9.297213e-04
## weight       -1.046022e-03
## cylinders    -3.826930e-01
## horsepower    7.100215e-05
pred.lda <- predict(autos.lda, test, type="response")
matriz <- confusionMatrix(pred.lda$class, as.factor(test$mpg01))
matriz$table
##           Reference
## Prediction  0  1
##          0 42  1
##          1  7 48
imprimirASS(matriz)
## [1] "Accuracy: 0.9184 Sensitivity: 0.8571 Specificity: 0.9796"

Podemos apreciar que la tasa de error es de un 10.2%

autos.qda <- qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = train)
autos.qda
## Call:
## qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = train)
## 
## Prior probabilities of groups:
##   0   1 
## 0.5 0.5 
## 
## Group means:
##   cylinders   weight displacement horsepower
## 0  6.789116 3601.537     271.4830  129.09524
## 1  4.217687 2352.027     118.0782   79.47619
pred.qda <- predict(autos.qda, test)
matriz <- confusionMatrix(pred.qda$class, as.factor(test$mpg01))
matriz$table
##           Reference
## Prediction  0  1
##          0 43  3
##          1  6 46
imprimirASS(matriz)
## [1] "Accuracy: 0.9082 Sensitivity: 0.8776 Specificity: 0.9388"

Como podemos ver la tasa de error es del 11.22%

autos.glm <- glm(mpg01 ~ cylinders + weight + displacement + horsepower, data = train, family = binomial)
summary(autos.glm)
## 
## Call:
## glm(formula = mpg01 ~ cylinders + weight + displacement + horsepower, 
##     family = binomial, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3670  -0.1738   0.0466   0.3810   3.1933  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  11.7659288  1.9131492   6.150 7.75e-10 ***
## cylinders    -0.0161152  0.3729110  -0.043  0.96553    
## weight       -0.0020281  0.0008052  -2.519  0.01177 *  
## displacement -0.0092265  0.0090618  -1.018  0.30859    
## horsepower   -0.0447260  0.0155085  -2.884  0.00393 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 407.57  on 293  degrees of freedom
## Residual deviance: 163.52  on 289  degrees of freedom
## AIC: 173.52
## 
## Number of Fisher Scoring iterations: 7
p <- predict(autos.glm, test, type = "response")
predichos=as.factor(ifelse(test = p > 0.5, yes = 1, no = 0))
verdaderos=as.factor(test$mpg01)
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 43  1
##          1  6 48
## [1] "Accuracy: 0.9286 Sensitivity: 0.8776 Specificity: 0.9796"

La tasa de error es de 11.22%

train.X2 = cbind(train$displacement, train$weight, train$cylinders, train$year)
test.X2 = cbind(test$displacement, test$weight, test$cylinders, test$year)
train.Y2 = cbind(train$mpg01)
# Para K = 10
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 10)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table
##           Reference
## Prediction  0  1
##          0 44  2
##          1  5 47
## [1] "Accuracy: 0.9286 Sensitivity: 0.898 Specificity: 0.9592"

La tasa de error el de 10.2%

# Para K = 50
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 50)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table
##           Reference
## Prediction  0  1
##          0 43  3
##          1  6 46
## [1] "Accuracy: 0.9082 Sensitivity: 0.8776 Specificity: 0.9388"

La tasa de error es del 10.2%

# Para K = 100
pred.knn <- knn(train.X2, test.X2, train.Y2, k = 100)
z <-confusionMatrix(pred.knn, as.factor(test$mpg01))
z$table
##           Reference
## Prediction  0  1
##          0 42  3
##          1  7 46
## [1] "Accuracy: 0.898 Sensitivity: 0.8571 Specificity: 0.9388"

La tasa de error es del 10.2%

Podemos ver que a partir de k = 10, la tasa de error se mantiene.

12

Power <- function() {
  2^3
}
Power()
## [1] 8
Power2 <- function(x, a) {
  x^a
}
Power2(3,8)
## [1] 6561
Power2(10,3)
## [1] 1000
Power2(8,17)
## [1] 2.2518e+15
Power2(131,3)
## [1] 2248091
Power3 <- function(x , a) {
    result <- x^a
    return(result)
}
x <- 1:10
plot(x, Power3(x, 2), log = "xy", xlab = "Log de x", ylab = "Log de x^2", main = "Log de x^2 vs Log de x")

PlotPower <- function(x, a) {
    plot(x, Power3(x, a))
}

PlotPower(1:10, 3)

13

Se requiere usar la base de datos “Boston”, la cual a continuación podemos apreciar las variables que posee:

str(Boston)
## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

Ahora calcularemos la mediana de la variiable crim para crear una nueva variable, la cual será binaria en donde 0 = crim <= mediana ó 1 = crim > mediana

mediana_crim <- rep(0, length(Boston$crim))
mediana_crim[Boston$crim > median(Boston$crim)] <- 1
Boston$mediana_crim = mediana_crim
head(Boston)
##      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12
##   lstat medv mediana_crim
## 1  4.98 24.0            0
## 2  9.14 21.6            0
## 3  4.03 34.7            0
## 4  2.94 33.4            0
## 5  5.33 36.2            0
## 6  5.21 28.7            0

Ahora construiremos nuestros datos de entrenamiento y de prueba

set.seed(123)
train_index <- sample(1:nrow(Boston), 0.8 * nrow(Boston))
test_index <- setdiff(1:nrow(Boston), train_index)
train <- Boston[train_index,0:length(Boston)]
test <- Boston[test_index,0:length(Boston)]

Ahora para entender un poco mejor las variables realizaremos un gráfico de dispersión y de correlaciones

pairs.panels(Boston, 
             method = "pearson",
             hist.col = "#00AFBB",
             density = TRUE,  
             )

Podemos apreciar que respecto a la variable mediana_crim, las variables nox, rad, dis, age, tax e indus son las que más están correlacionadas. Para este ejercicio se utilizarán modelos partiendo de una variable (la de mayor correlación) y se irán agregando nuevas variables para determinar cual es el mejor modelo para cada ténica.

Comenzando con modelos de regresión logística:

  • Modelo con nox como única variable predictora
glm_13.fit <- glm(mediana_crim ~ nox, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
verdaderos=as.factor(test[,15])
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 40  5
##          1  9 48
## [1] "Accuracy: 0.8627 Sensitivity: 0.8163 Specificity: 0.9057"
  • Modelo con nox y rad como únicas variables predictoras
glm_13.fit <- glm(mediana_crim ~ nox + rad, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 46  7
##          1  3 46
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
  • Modelo con nox, rad y dis como únicas variables predictoras
glm_13.fit <- glm(mediana_crim ~ nox + rad + dis, data = train, family = binomial)
probs <- predict(glm_13.fit, test[,-15], type = "response")
predichos=as.factor(ifelse(test = probs > 0.5, yes = 1, no = 0))
matriz<-confusionMatrix(predichos, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 46  7
##          1  3 46
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"

Como podemos apreciar en los anteriores modelos, con las variables nox y rad como únicas predictoras se obtuvo el mejor modelo, cuando agregamos más variables predictoras no se mejoran los resultados: Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679

Ahora seguimos con los modelos LDA

  • Modelo con nox como única variable predictora
ldapre <- lda(mediana_crim ~ nox, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 42 11
##          1  7 42
## [1] "Accuracy: 0.8235 Sensitivity: 0.8571 Specificity: 0.7925"
  • Modelo con nox y rad como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 47 12
##          1  2 41
## [1] "Accuracy: 0.8627 Sensitivity: 0.9592 Specificity: 0.7736"
  • Modelo con nox, rad y dis como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad + dis, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 45 11
##          1  4 42
## [1] "Accuracy: 0.8529 Sensitivity: 0.9184 Specificity: 0.7925"
  • Modelo con nox, rad, dis y age como únicas variables predictoras
ldapre <- lda(mediana_crim ~ nox + rad + dis + age, data = train)
pred.lda <- predict(ldapre, test[,-15])
matriz<-confusionMatrix(pred.lda$class, verdaderos)
matriz$table
##           Reference
## Prediction  0  1
##          0 45 13
##          1  4 40
## [1] "Accuracy: 0.8333 Sensitivity: 0.9184 Specificity: 0.7547"

Podemos ver que aumentando el número de variables no presenta un mejora en la precisión de nuestro modelo, de lo que se probaron el mejor fue el que se usó con las variables nox y rad con los siguientes parámetros: Accuracy: 0.8627 Sensitivity: 0.9592 Specificity: 0.7736

Modelos KNN

modelos_knn <- function(x_train, x_test, y_train, y_test, nk){
  pred.knn <- knn(x_train, x_test, y_train, k = nk)
  z <-confusionMatrix(pred.knn, y_test)
  print(paste0("K = ", nk))
  imprimirASS(z)
}
entrenamiento.x <- as.matrix(train[,-15])
prueba.x <- as.matrix(test[,-15])
entrenamiento.y <- as.factor(train[,15])
prueba.y <- as.factor(test[,15])
for (k in c(2:20)) {
  modelos_knn(entrenamiento.x, prueba.x, entrenamiento.y, prueba.y, k)
}
## [1] "K = 2"
## [1] "Accuracy: 0.9118 Sensitivity: 0.9184 Specificity: 0.9057"
## [1] "K = 3"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9184 Specificity: 0.9245"
## [1] "K = 4"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9388 Specificity: 0.9245"
## [1] "K = 5"
## [1] "Accuracy: 0.951 Sensitivity: 0.9796 Specificity: 0.9245"
## [1] "K = 6"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 7"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 8"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 9"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 10"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 11"
## [1] "Accuracy: 0.9314 Sensitivity: 0.9592 Specificity: 0.9057"
## [1] "K = 12"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9592 Specificity: 0.8868"
## [1] "K = 13"
## [1] "Accuracy: 0.9118 Sensitivity: 0.9388 Specificity: 0.8868"
## [1] "K = 14"
## [1] "Accuracy: 0.9216 Sensitivity: 0.9388 Specificity: 0.9057"
## [1] "K = 15"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 16"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 17"
## [1] "Accuracy: 0.902 Sensitivity: 0.9388 Specificity: 0.8679"
## [1] "K = 18"
## [1] "Accuracy: 0.8922 Sensitivity: 0.9184 Specificity: 0.8679"
## [1] "K = 19"
## [1] "Accuracy: 0.8725 Sensitivity: 0.9184 Specificity: 0.8302"
## [1] "K = 20"
## [1] "Accuracy: 0.8824 Sensitivity: 0.9184 Specificity: 0.8491"

Podemos ver con un K = 5 se obtuvo la mayor precisión, lo cual también determina, que éste es el mejor modelo de todos los que se probaron en el ejercicio

Capitulo 8

7

library(MASS)
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:psych':
## 
##     outlier
set.seed(1)

train <- sample(1:nrow(Boston), 400)
Boston.train <- Boston[train, -14]
Boston.test <- Boston[-train, -14]
Y.train <- Boston[train, 14]
Y.test <- Boston[-train, 14]
boston1 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = ncol(Boston) - 1, ntree = 500)
boston2 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = (ncol(Boston) - 1) / 2, ntree = 500)
boston3 <- randomForest(Boston.train, y = Y.train, xtest = Boston.test, ytest = Y.test, mtry = sqrt(ncol(Boston) - 1), ntree = 500)

valores <- data.frame(c(1:500))
names(valores) <- "arboles"
valores["MSE1"] <- boston1$test$mse
valores["MSE2"] <- boston2$test$mse
valores["MSE3"] <- boston3$test$mse

library(ggplot2)
library(reshape2)
dd = melt(valores, id=c("arboles"))

theme_set(theme_bw()) 
ggplot(dd) + geom_line(aes(x=arboles, y=value, color=variable)) +
  scale_color_manual(name = "m = ",
                     labels = c("p","p/2","\u221Ap"),
                     values=c("green","red","blue")) +
  labs(x = "Número de árboles",
       y = "Error de clasificación") +
  theme(legend.position="top")

Como podemos apreciar en el gráfico anterior con solamente un árbol el erro tiende a ser muy alto, pero a medida que se aumenta el número de árboles, se estabiliza este error, aproximadamente a partir de 100 árboles, este comportamiento se presenta en general con los tres valores de m. Respecto a con cual m se presenta el menor error, podemos ver que en mi caso fue con \(m=p\), aunque si cambiamos la semilla este resultado puede cambiar.

8

train <- sample(1:nrow(Carseats), 300)
Carseats.train <- Carseats[train, ]
Carseats.test <- Carseats[-train, ]
library(rpart)
library(rpart.plot)

tree <- rpart(Sales~., data=Carseats.train)

rpart.plot(tree, box.palette="RdBu", nn=TRUE)

yhat <- predict(tree, newdata = Carseats.test)
m1 <- mean((yhat - Carseats.test$Sales)^2)
m1
## [1] 4.197228

Como podemos apreciar el error MSE de la prueba es de 4.197228

printcp(tree)
## 
## Regression tree:
## rpart(formula = Sales ~ ., data = Carseats.train)
## 
## Variables actually used in tree construction:
## [1] Advertising Age         CompPrice   Income      Population  Price      
## [7] ShelveLoc  
## 
## Root node error: 2247.1/300 = 7.4904
## 
## n= 300 
## 
##          CP nsplit rel error  xerror     xstd
## 1  0.211723      0   1.00000 1.01804 0.082943
## 2  0.123518      1   0.78828 0.80759 0.065671
## 3  0.057792      2   0.66476 0.71205 0.057322
## 4  0.039590      3   0.60697 0.68368 0.053027
## 5  0.032717      4   0.56738 0.67813 0.053775
## 6  0.026549      5   0.53466 0.66341 0.050763
## 7  0.021459      6   0.50811 0.62629 0.048310
## 8  0.018375      7   0.48665 0.63310 0.050340
## 9  0.017637      9   0.44990 0.64361 0.051830
## 10 0.017500     10   0.43226 0.64903 0.051778
## 11 0.015287     11   0.41476 0.62861 0.050278
## 12 0.014037     13   0.38419 0.62867 0.050364
## 13 0.013272     14   0.37015 0.64636 0.051480
## 14 0.012158     15   0.35688 0.64895 0.051655
## 15 0.011647     16   0.34472 0.64894 0.051919
## 16 0.010000     17   0.33308 0.64854 0.051978
plotcp(tree)

min <- which.min(tree$cptable[,"xerror"])
min
## 7 
## 7
cp <-tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]
cp
## [1] 0.02145892

Con la validación cruzada obtenemos un tamaño de arbol igual a 7 y un cp de 0.02145892

ptree<- prune(tree, cp= cp)
rpart.plot(ptree, box.palette="RdBu", nn=TRUE)

yhat <- predict(ptree, newdata = Carseats.test)
m2 <- mean((yhat - Carseats.test$Sales)^2)
m2
## [1] 4.54531

Podemos ver que se aumenta el MSE en nuestro caso después de la poda

bag.carseats <- randomForest(Sales ~ ., data = Carseats.train, mtry = (ncol(Carseats) - 1), importance = TRUE)
yhat.bag <- predict(bag.carseats, newdata = Carseats.test)
mean((yhat.bag - Carseats.test$Sales)^2)
## [1] 2.381242

Podemos ver como se mejoró a 2.38 el MSE con este enfoque.

importance(bag.carseats)
##               %IncMSE IncNodePurity
## CompPrice   38.928557    272.987960
## Income      11.725449    122.545479
## Advertising 14.937292    128.773274
## Population   2.518008     95.322579
## Price       71.644649    716.074862
## ShelveLoc   69.101715    596.025736
## Age         17.537409    167.453663
## Education    2.793162     64.343849
## Urban       -2.073262      8.985155
## US           4.206136      9.651551

Y con esto podemos apreciar que Price, ShelveLoc y CompPrice son las variables predictoras más importantes en este análisis

bag.carseats <- randomForest(Sales ~ ., data = Carseats.train, mtry = sqrt((ncol(Carseats) - 1)), importance = TRUE)
yhat.bag <- predict(bag.carseats, newdata = Carseats.test)
mean((yhat.bag - Carseats.test$Sales)^2)
## [1] 2.80204

En este caso tenemos un MSE de 2.8

importance(bag.carseats)
##                %IncMSE IncNodePurity
## CompPrice   16.2755940     216.47915
## Income       7.8217396     184.80927
## Advertising 14.1807157     183.03683
## Population   0.6892522     155.21205
## Price       42.0256227     551.81898
## ShelveLoc   47.8554474     473.99364
## Age         12.1305637     213.32128
## Education    2.8662179      91.35868
## Urban       -2.3838218      18.20191
## US           3.4288791      26.69131

También podemos ver que las variables más importantes siguen siendo las mismas que las que se declararon en el punto anterior.

9

train <- sample(1:nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test <- OJ[-train, ]
library(tree)
tree.oj <- tree(Purchase ~ ., data = OJ.train)
summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH"     "SalePriceMM" "SpecialCH"   "PriceDiff"  
## Number of terminal nodes:  9 
## Residual mean deviance:  0.707 = 559.3 / 791 
## Misclassification error rate: 0.1512 = 121 / 800

Como podemos apreciar, el arbol tiene 9 nodos terminales y un error de entrenamiento de 0.1512. También podemos ver que las variables usadas en la construcción del árbol son LoyalCH, SalePriceMM, specialCH y PriceDiff

tree.oj
## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 800 1065.00 CH ( 0.61625 0.38375 )  
##    2) LoyalCH < 0.5036 351  424.90 MM ( 0.29345 0.70655 )  
##      4) LoyalCH < 0.264232 151  102.10 MM ( 0.10596 0.89404 )  
##        8) LoyalCH < 0.051325 57    0.00 MM ( 0.00000 1.00000 ) *
##        9) LoyalCH > 0.051325 94   85.77 MM ( 0.17021 0.82979 ) *
##      5) LoyalCH > 0.264232 200  273.90 MM ( 0.43500 0.56500 )  
##       10) SalePriceMM < 2.04 109  126.30 MM ( 0.26606 0.73394 )  
##         20) SpecialCH < 0.5 83   78.43 MM ( 0.18072 0.81928 ) *
##         21) SpecialCH > 0.5 26   35.89 CH ( 0.53846 0.46154 ) *
##       11) SalePriceMM > 2.04 91  119.20 CH ( 0.63736 0.36264 ) *
##    3) LoyalCH > 0.5036 449  349.40 CH ( 0.86860 0.13140 )  
##      6) LoyalCH < 0.764572 188  217.80 CH ( 0.73404 0.26596 )  
##       12) PriceDiff < 0.265 113  153.40 CH ( 0.58407 0.41593 )  
##         24) PriceDiff < -0.165 34   41.19 MM ( 0.29412 0.70588 ) *
##         25) PriceDiff > -0.165 79   95.30 CH ( 0.70886 0.29114 ) *
##       13) PriceDiff > 0.265 75   25.19 CH ( 0.96000 0.04000 ) *
##      7) LoyalCH > 0.764572 261   78.30 CH ( 0.96552 0.03448 ) *

Se escoge el nodo terminal 7, en este podemos apreciar que su creterio de división es LoyalCH > 0.764572, éste tiene 261 observaciones en este nodo con una desviación de 78.30. La predicción en este nodo es Sales = CH, también podemos ver que el 96.552% de las observaciones toman el valor de CH, mientras que el 3.448% toman el valor de MM.

plot(tree.oj)
text(tree.oj, pretty = 0)

En este caso podemos apreciar que la variable principal es LoyalCH, los nodos que nacen de este también usan LoyalCH para dividir el arbol. Podemos ver que en el nodo principal si LoyalCH es mayor o igual a 0.5036, sólo se usa PriceDiff como otra variable, mientras en caso contrario se utilizan las variables salePriceMM y SpecialCH.

tree.pred <- predict(tree.oj, OJ.test, type = "class")

matriz<-confusionMatrix(tree.pred, OJ.test$Purchase)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  CH  MM
##         CH 144  34
##         MM  16  76
##                                           
##                Accuracy : 0.8148          
##                  95% CI : (0.7633, 0.8593)
##     No Information Rate : 0.5926          
##     P-Value [Acc > NIR] : 4.577e-15       
##                                           
##                   Kappa : 0.6064          
##                                           
##  Mcnemar's Test P-Value : 0.01621         
##                                           
##             Sensitivity : 0.9000          
##             Specificity : 0.6909          
##          Pos Pred Value : 0.8090          
##          Neg Pred Value : 0.8261          
##              Prevalence : 0.5926          
##          Detection Rate : 0.5333          
##    Detection Prevalence : 0.6593          
##       Balanced Accuracy : 0.7955          
##                                           
##        'Positive' Class : CH              
## 

Podemos ver que la precisión del modelo es de 81.48%, siendo mucho más preciso para predecir CH con un 90% mientras que para predecir MM sólo es de un 69.09%.

cv.oj <- cv.tree(tree.oj, FUN = prune.misclass)
cv.oj
## $size
## [1] 9 8 7 4 2 1
## 
## $dev
## [1] 135 135 140 151 172 303
## 
## $k
## [1]       -Inf   0.000000   2.000000   4.666667  12.500000 145.000000
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"
valores <- data.frame(cv.oj$size, cv.oj$dev)
names(valores) <- c("Size", "Dev")

library(ggplot2)
library(reshape2)

theme_set(theme_bw()) 
ggplot(data = valores, aes(x=Size, y=Dev)) + 
  geom_line(color="red", linetype = "dashed") +
  geom_point() +
  scale_x_discrete(limits=c(1:9))

El tamaño de árbol con menor error es de 7

prune.oj <- prune.misclass(tree.oj, best = 7)
plot(prune.oj)
text(prune.oj, pretty = 0)

summary(prune.oj)
## 
## Classification tree:
## snip.tree(tree = tree.oj, nodes = c(4L, 10L))
## Variables actually used in tree construction:
## [1] "LoyalCH"     "SalePriceMM" "PriceDiff"  
## Number of terminal nodes:  7 
## Residual mean deviance:  0.7409 = 587.5 / 793 
## Misclassification error rate: 0.1538 = 123 / 800
summary(tree.oj)
## 
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH"     "SalePriceMM" "SpecialCH"   "PriceDiff"  
## Number of terminal nodes:  9 
## Residual mean deviance:  0.707 = 559.3 / 791 
## Misclassification error rate: 0.1512 = 121 / 800

Vemos que practicamente no existe diferencia entre los dos árboles, siendo levemenete mayor en el arbol podado

prune.pred <- predict(prune.oj, OJ.test, type = "class")

matriz<-confusionMatrix(prune.pred, OJ.test$Purchase)
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  CH  MM
##         CH 141  34
##         MM  19  76
##                                           
##                Accuracy : 0.8037          
##                  95% CI : (0.7512, 0.8494)
##     No Information Rate : 0.5926          
##     P-Value [Acc > NIR] : 1.152e-13       
##                                           
##                   Kappa : 0.5846          
##                                           
##  Mcnemar's Test P-Value : 0.05447         
##                                           
##             Sensitivity : 0.8812          
##             Specificity : 0.6909          
##          Pos Pred Value : 0.8057          
##          Neg Pred Value : 0.8000          
##              Prevalence : 0.5926          
##          Detection Rate : 0.5222          
##    Detection Prevalence : 0.6481          
##       Balanced Accuracy : 0.7861          
##                                           
##        'Positive' Class : CH              
## 

Vemos que bajó la precisión del modelo con respecto al árbol sin podar.

10

Hitters <- na.omit(Hitters)
Hitters$Salary <- log(Hitters$Salary)
train <- 1:200
Hitters.train = Hitters[train, ]
Hitters.test = Hitters[-train, ]
library(gbm)
## Loaded gbm 2.1.5
lambdas = seq(0.001, 0.3, 0.005)
mse <- rep(0, length(lambdas))
for (i in 1:length(lambdas)) {
  boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
    pred.train <- predict(boost.hitters, Hitters.train, n.trees = 1000)
    mse[i] <- mean((pred.train - Hitters.train$Salary)^2)
}
library(ggplot2)

mse_graph <- data.frame(lambdas)
mse_graph["MSE"]<- mse
names(mse_graph) <- c("Lambdas", "MSE")

theme_set(theme_bw()) 
ggplot(data = mse_graph, aes(x = Lambdas, y = MSE)) +
  geom_line(color="red", linetype = "dashed") +
  geom_point() 

mse2 <- rep(NA, length(lambdas))
for (i in 1:length(lambdas)) {
    boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
    yhat <- predict(boost.hitters, Hitters.test, n.trees = 1000)
    mse2[i] <- mean((yhat - Hitters.test$Salary)^2)
}
library(ggplot2)

mse_graph2 <- data.frame(lambdas)
mse_graph2["MSE"]<- mse2
names(mse_graph2) <- c("Lambdas", "MSE")

theme_set(theme_bw()) 
ggplot(data = mse_graph2, aes(x = Lambdas, y = MSE)) +
  geom_line(color="red", linetype = "dashed") +
  geom_point() 

min(mse2)
## [1] 0.2475835
lambdas[which.min(mse2)]
## [1] 0.166
lm810 <- lm(Salary ~ ., data = Hitters.train)
pred810 <- predict(lm810, Hitters.test)
mean((pred810 - Hitters.test$Salary)^2)
## [1] 0.4917959
library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
## 
##     R2
## The following object is masked from 'package:stats':
## 
##     loadings
pcr.fit=pcr(Salary~., data=Hitters.train , scale=TRUE , validation ="CV")
validationplot(pcr.fit ,val.type="MSEP")

pcr.pred=predict (pcr.fit ,Hitters.test,ncomp =1)
mean((pcr.pred -Hitters.test$Salary)^2)
## [1] 0.4661183

Podemos ver que el MSE por boosting es menor que el dado por la regresión lineal y componentes principales.

boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[which.min(mse2)])
summary(boost.hitters)

##                 var    rel.inf
## CAtBat       CAtBat 17.7715143
## CWalks       CWalks  9.0983419
## CRuns         CRuns  8.1091446
## PutOuts     PutOuts  7.8019290
## Walks         Walks  7.3620230
## CHmRun       CHmRun  6.4505765
## Assists     Assists  5.6779659
## Hits           Hits  5.0562630
## RBI             RBI  4.9022635
## Years         Years  4.8919157
## AtBat         AtBat  4.3145631
## CHits         CHits  4.0969282
## CRBI           CRBI  3.8001772
## Runs           Runs  3.6363269
## HmRun         HmRun  3.3847743
## Errors       Errors  2.3407756
## NewLeague NewLeague  0.5117111
## Division   Division  0.4565112
## League       League  0.3362948

La variable CatBat es la más importante.

bag.hitters <- randomForest(Salary ~ ., data = Hitters.train, mtry = 19)
yhat.bag <- predict(bag.hitters, newdata = Hitters.test)
mean((yhat.bag - Hitters.test$Salary)^2)
## [1] 0.22971

Podemos ver que el MSE fue de 0.007 lo cual significa que es mucho mejor que los probados hasta ahora.

11

train <- 1:1000
Caravan$Purchase <- ifelse(Caravan$Purchase == "Yes", 1, 0)
Caravan.train <- Caravan[train, ]
Caravan.test <- Caravan[-train, ]
bc <- gbm(Purchase ~ ., data = Caravan.train, distribution = "gaussian", n.trees = 1000, shrinkage = 0.01)
## Warning in gbm.fit(x = x, y = y, offset = offset, distribution =
## distribution, : variable 50: PVRAAUT has no variation.
## Warning in gbm.fit(x = x, y = y, offset = offset, distribution =
## distribution, : variable 71: AVRAAUT has no variation.
summary(bc)

##               var    rel.inf
## PPERSAUT PPERSAUT 14.6627387
## MKOOPKLA MKOOPKLA 11.5135457
## MOPLHOOG MOPLHOOG  8.1166824
## MBERMIDD MBERMIDD  5.5659055
## PBRAND     PBRAND  4.8848882
## MGODGE     MGODGE  4.3094073
## ABRAND     ABRAND  4.0899213
## MINK3045 MINK3045  3.7799757
## MAUT1       MAUT1  3.0595852
## MOSTYPE   MOSTYPE  2.6279165
## PWAPART   PWAPART  2.6176353
## MBERARBG MBERARBG  2.1522326
## MGODPR     MGODPR  2.1053373
## MAUT2       MAUT2  2.0837180
## MSKC         MSKC  2.0548149
## MSKA         MSKA  1.8052960
## MSKB1       MSKB1  1.7245921
## PBYSTAND PBYSTAND  1.5426299
## MBERHOOG MBERHOOG  1.4526150
## MFWEKIND MFWEKIND  1.4470476
## MINK7512 MINK7512  1.3601434
## MGODRK     MGODRK  1.3247916
## MGODOV     MGODOV  1.2359721
## MOPLMIDD MOPLMIDD  1.1189140
## MINK4575 MINK4575  1.0021078
## MRELGE     MRELGE  0.9752463
## MRELOV     MRELOV  0.9168410
## MAUT0       MAUT0  0.8624586
## MINKGEM   MINKGEM  0.8250010
## MFGEKIND MFGEKIND  0.7750144
## APERSAUT APERSAUT  0.7452557
## MBERBOER MBERBOER  0.6269871
## MINKM30   MINKM30  0.6230703
## MOSHOOFD MOSHOOFD  0.6128366
## MZFONDS   MZFONDS  0.5781958
## PMOTSCO   PMOTSCO  0.5172908
## MHHUUR     MHHUUR  0.4691845
## MHKOOP     MHKOOP  0.4467587
## MZPART     MZPART  0.4351830
## MBERARBO MBERARBO  0.4270172
## MGEMLEEF MGEMLEEF  0.3961990
## MGEMOMV   MGEMOMV  0.3886987
## MSKD         MSKD  0.3409863
## PLEVEN     PLEVEN  0.3215526
## MSKB2       MSKB2  0.2531789
## MINK123M MINK123M  0.2416957
## MOPLLAAG MOPLLAAG  0.1897649
## MRELSA     MRELSA  0.1837443
## MBERZELF MBERZELF  0.1086964
## MFALLEEN MFALLEEN  0.1007280
## MAANTHUI MAANTHUI  0.0000000
## PWABEDR   PWABEDR  0.0000000
## PWALAND   PWALAND  0.0000000
## PBESAUT   PBESAUT  0.0000000
## PVRAAUT   PVRAAUT  0.0000000
## PAANHANG PAANHANG  0.0000000
## PTRACTOR PTRACTOR  0.0000000
## PWERKT     PWERKT  0.0000000
## PBROM       PBROM  0.0000000
## PPERSONG PPERSONG  0.0000000
## PGEZONG   PGEZONG  0.0000000
## PWAOREG   PWAOREG  0.0000000
## PZEILPL   PZEILPL  0.0000000
## PPLEZIER PPLEZIER  0.0000000
## PFIETS     PFIETS  0.0000000
## PINBOED   PINBOED  0.0000000
## AWAPART   AWAPART  0.0000000
## AWABEDR   AWABEDR  0.0000000
## AWALAND   AWALAND  0.0000000
## ABESAUT   ABESAUT  0.0000000
## AMOTSCO   AMOTSCO  0.0000000
## AVRAAUT   AVRAAUT  0.0000000
## AAANHANG AAANHANG  0.0000000
## ATRACTOR ATRACTOR  0.0000000
## AWERKT     AWERKT  0.0000000
## ABROM       ABROM  0.0000000
## ALEVEN     ALEVEN  0.0000000
## APERSONG APERSONG  0.0000000
## AGEZONG   AGEZONG  0.0000000
## AWAOREG   AWAOREG  0.0000000
## AZEILPL   AZEILPL  0.0000000
## APLEZIER APLEZIER  0.0000000
## AFIETS     AFIETS  0.0000000
## AINBOED   AINBOED  0.0000000
## ABYSTAND ABYSTAND  0.0000000

Las variables PPERSAUT y MKOOPKLA son las más importantes

probs.test <- predict(bc, Caravan.test, n.trees = 1000, type = "response")
pred.test <- ifelse(probs.test > 0.2, 1, 0)
matriz <- confusionMatrix(as.factor(pred.test), as.factor(Caravan.test$Purchase))
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 4495  277
##          1   38   12
##                                           
##                Accuracy : 0.9347          
##                  95% CI : (0.9273, 0.9415)
##     No Information Rate : 0.9401          
##     P-Value [Acc > NIR] : 0.9446          
##                                           
##                   Kappa : 0.0541          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.99162         
##             Specificity : 0.04152         
##          Pos Pred Value : 0.94195         
##          Neg Pred Value : 0.24000         
##              Prevalence : 0.94007         
##          Detection Rate : 0.93219         
##    Detection Prevalence : 0.98963         
##       Balanced Accuracy : 0.51657         
##                                           
##        'Positive' Class : 0               
## 

Vemos que la predicción de personas que en realidad harán la compra tiene una precisión del 3.46%

logit.caravan <- glm(Purchase ~ ., data = Caravan.train, family = "binomial")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
probs.test2 <- predict(logit.caravan, Caravan.test, type = "response")
## Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type
## == : prediction from a rank-deficient fit may be misleading
pred.test2 <- ifelse(probs.test > 0.2, 1, 0)
matriz <- confusionMatrix(as.factor(pred.test2), as.factor(Caravan.test$Purchase))
matriz
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 4495  277
##          1   38   12
##                                           
##                Accuracy : 0.9347          
##                  95% CI : (0.9273, 0.9415)
##     No Information Rate : 0.9401          
##     P-Value [Acc > NIR] : 0.9446          
##                                           
##                   Kappa : 0.0541          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.99162         
##             Specificity : 0.04152         
##          Pos Pred Value : 0.94195         
##          Neg Pred Value : 0.24000         
##              Prevalence : 0.94007         
##          Detection Rate : 0.93219         
##    Detection Prevalence : 0.98963         
##       Balanced Accuracy : 0.51657         
##                                           
##        'Positive' Class : 0               
## 

Vemos que usando una regresión logística la precisión del modelo a la hora de predecir personas que realmente realizaron la compra es de un 3.46%

12

Se usará la base de datos ‘Smarket’, el cual representa los porcentajes diarios de rendimiento para el índice bursátil S&P 500 entre 2001 y 2005. En este se intentará predecir la variable ‘Direction’

train <- sample(nrow(Weekly), nrow(Weekly)*0.7)
Weekly$Direction <- ifelse(Weekly$Direction == "Up", 1, 0)
Weekly.train <- Weekly[train, ]
Weekly.test <- Weekly[-train, ]
  • Aplicando boosting
bf <- gbm(Direction ~ . - Year - Today, data = Weekly.train, distribution = "bernoulli", n.trees = 5000)
bprobs <- predict(bf, newdata = Weekly.test, n.trees = 5000)
bpred <- ifelse(bprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(bpred), as.factor(Weekly.test$Direction))
matriz$table
##           Reference
## Prediction   0   1
##          0  84 100
##          1  54  89
## [1] "Accuracy: 0.5291 Sensitivity: 0.6087 Specificity: 0.4709"
  • Aplicando bagging
bagf <- randomForest(Direction ~ . - Year - Today, data = Weekly.train, mtry = 6)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?
bagprobs <- predict(bagf, newdata = Weekly.test)
bagpred <- ifelse(bagprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(bagpred), as.factor(Weekly.test$Direction))
matriz$table
##           Reference
## Prediction   0   1
##          0  56  57
##          1  82 132
## [1] "Accuracy: 0.5749 Sensitivity: 0.4058 Specificity: 0.6984"
  • Aplicando random Forest
rff <- randomForest(Direction ~ . - Year - Today, data = Weekly.train, mtry = 2)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values. Are you sure you want to do regression?
rfprobs <- predict(rff, newdata = Weekly.test)
rfpred <- ifelse(rfprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(rfpred), as.factor(Weekly.test$Direction))
matriz$table
##           Reference
## Prediction   0   1
##          0  58  54
##          1  80 135
## [1] "Accuracy: 0.5902 Sensitivity: 0.4203 Specificity: 0.7143"
  • Regresión logística
logf <- glm(Direction ~ . - Year - Today, data = Weekly.train, family = "binomial")
logprobs <- predict(logf, newdata = Weekly.test, type = "response")
logpred <- ifelse(logprobs > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(logpred), as.factor(Weekly.test$Direction))
matriz$table
##           Reference
## Prediction   0   1
##          0   7   8
##          1 131 181
## [1] "Accuracy: 0.5749 Sensitivity: 0.0507 Specificity: 0.9577"
  • Regresión lineal
lmf <- lm(Direction ~ . - Year - Today, data = Weekly.train)
lmprob <- predict(lmf, Weekly.test)
lmpred <- ifelse(lmprob > 0.5, 1, 0)
matriz <- confusionMatrix(as.factor(lmpred), as.factor(Weekly.test$Direction))
matriz$table
##           Reference
## Prediction   0   1
##          0   7   8
##          1 131 181
## [1] "Accuracy: 0.5749 Sensitivity: 0.0507 Specificity: 0.9577"

En general no se encuentran grandes diferencias entre ellos, pero de estos propuestos, el que mejor resultados da es el bagging con un accuraccy de 54.74%

Capitulo 9

4

df <- data.frame(replicate(2, rnorm(n = 100)))
df <- as.tibble(df)
## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.
gen <- function(x, y) {
    x^2 + y^2 <= 1
}

df <- df %>%
    rename(Var1 = X1, Var2 = X2) %>%
    mutate(Clases = ifelse(gen(Var1, Var2),
                          'Clase 1',
                          'Clase 2'),
           Clases = factor(Clases))

inTrain <- sample(nrow(df), nrow(df)*0.7)
train <- df[inTrain,]
test <- df[-inTrain,]

theme_set(theme_bw())
ggplot(df, aes(Var1, Var2, col = Clases)) +
    geom_point(size = 2)

  • Clasificador de soporte vectorial
svm1 <- svm(Clases ~ ., data = train, 
                 kernel = 'linear', 
                 scale = FALSE, cost = 10)
plot(svm1, data = train)

matriz <- confusionMatrix(predict(svm1, test), test$Clases)
matriz$table
##           Reference
## Prediction Clase 1 Clase 2
##    Clase 1       8      10
##    Clase 2       4       8
## [1] "Accuracy: 0.5333 Sensitivity: 0.6667 Specificity: 0.4444"
plot(svm1, data = test)

  • SVM con kernel polinomial
svm2 <- svm(Clases ~ ., data = train, 
                kernel = 'polynomial', degree = 2,
                scale = FALSE, cost = 1)

plot(svm2, train)

matriz <- confusionMatrix(predict(svm2, test), test$Clases)
matriz$table
##           Reference
## Prediction Clase 1 Clase 2
##    Clase 1      12       1
##    Clase 2       0      17
## [1] "Accuracy: 0.9667 Sensitivity: 1 Specificity: 0.9444"
plot(svm2, data = test)

  • SVM con kernel radial
svm3 <- svm(Clases ~ ., data = train,
                  kernel = 'radial',
                  scale = FALSE, cost = 1)
plot(svm3, train)

matriz <- confusionMatrix(predict(svm3, test), test$Clases)
matriz$table
##           Reference
## Prediction Clase 1 Clase 2
##    Clase 1      11       2
##    Clase 2       1      16
## [1] "Accuracy: 0.9 Sensitivity: 0.9167 Specificity: 0.8889"
plot(svm3, data = test)

Podemos apreciar que, tal como lo decía el enunciado, tanto la SVM de kernel radial como polinomial tienen un mejor rendimiento que el modelo con kernel lineal. El modelo que tuvo un mejor rendimiento en este experimento fue la SVM polinomial de grado 2, con una precisión del 100% con el set de entrenamiento.

5

x1 <- runif(500) - 0.5
x2 <- runif(500) - 0.5
y <- 1 * (x1^2 - x2^2 > 0)
df <- data.frame(y,x1,x2)

df <- df %>%
    mutate(y = factor(y ))

theme_set(theme_bw()) 
ggplot(df, aes(x=x1, y=x2, col=y)) + 
  geom_point()

train <- df
logreg_fit <- glm(y ~ ., data = train, family = 'binomial')
summary(logreg_fit)
## 
## Call:
## glm(formula = y ~ ., family = "binomial", data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3452  -1.0996  -0.9886   1.1824   1.3738  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -0.038121   0.090234  -0.422  0.67269   
## x1           0.008842   0.308455   0.029  0.97713   
## x2           0.857039   0.309664   2.768  0.00565 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 692.86  on 499  degrees of freedom
## Residual deviance: 685.07  on 497  degrees of freedom
## AIC: 691.07
## 
## Number of Fisher Scoring iterations: 4
probs <- predict(logreg_fit, data = train, type = "response")
preds <- rep(0, 500)
preds[probs > 0.5] <- 1
train['preds'] <- preds
theme_set(theme_bw()) 
train <- train %>%
    mutate(preds = factor(preds))
ggplot(train, aes(x=x1, y=x2, col=preds)) + 
  geom_point()

lm.fit = glm(y ~ poly(x1, 2) + x2, data = train, family = binomial)
summary(lm.fit)
## 
## Call:
## glm(formula = y ~ poly(x1, 2) + x2, family = binomial, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2464  -0.7807  -0.5253   0.7573   1.8684  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.05209    0.11049   0.471    0.637    
## poly(x1, 2)1  0.16859    2.62143   0.064    0.949    
## poly(x1, 2)2 32.04520    3.01261  10.637   <2e-16 ***
## x2            0.88434    0.36713   2.409    0.016 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 692.86  on 499  degrees of freedom
## Residual deviance: 516.82  on 496  degrees of freedom
## AIC: 524.82
## 
## Number of Fisher Scoring iterations: 4
probs <- predict(lm.fit, data = train, type = "response")
preds2 <- rep(0, 500)
preds2[probs > 0.5] <- 1
theme_set(theme_bw()) 
train <- train %>%
    mutate(preds2 = factor(preds2))
ggplot(train, aes(x=x1, y=x2, col=preds2)) + 
  geom_point()

svm.fit <- svm(y ~ x1 + x2, data = train, kernel = "linear")
preds3 <- predict(svm.fit, train)
theme_set(theme_bw()) 
train['preds3'] <- preds3 
ggplot(train, aes(x=x1, y=x2, col=preds3)) + 
  geom_point()

svmr <- svm(y ~ x1 + x2, data = train,
                  kernel = 'radial',
                  scale = FALSE)
predsr <- predict(svmr, train)
theme_set(theme_bw()) 
train['predsr'] <- predsr 
ggplot(train, aes(x=x1, y=x2, col=predsr)) + 
  geom_point()

Podemos ver a través de estos experimentos que las mquinas de soporte vectorial con metodos no lineas en datos que no tienen “perimetro” lineal, son muy efectivas, también si tienen perímetro lineal las svm con método lineal también son muy efectivas en estos casos.

6

x1=runif (1000) -0.5 
x2=runif (1000) -0.5

temp <- x1-x2

y <- rep(0, 1000)

for (i in 1:1000) {
  if (temp[i] > 0.1) {
    y[i] <- 1
  }
  else if (temp[i] < -0.1) {
    y[i] <- 0
  }
  else {
    y[i] <- runif(1) > 0.5
  }
}
  
df2 <- data.frame(x1,x2,y)

df2$y <- as.factor(df2$y)

ggplot(df2, aes(x=x1,y=x2,col=y)) +
  geom_point()

svm1 <- tune(svm, y ~ ., data = df2, kernel = "linear", ranges = list(cost = c(0.0001, 0.001, 0.01, 1, 5, 10, 100, 1000, 10000)))
summary(svm1)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##   cost
##  0.001
## 
## - best performance: 0.086 
## 
## - Detailed performance results:
##    cost error dispersion
## 1 1e-04 0.480 0.03915780
## 2 1e-03 0.086 0.04427189
## 3 1e-02 0.092 0.03910101
## 4 1e+00 0.089 0.03984693
## 5 5e+00 0.092 0.04417138
## 6 1e+01 0.092 0.04237400
## 7 1e+02 0.089 0.03813718
## 8 1e+03 0.089 0.03813718
## 9 1e+04 0.091 0.04040077

Podemos ver que un costo de 0.01 es el que nos dá el menor error

data.frame(cost = svm1$performance$cost, misclass = svm1$performance$error * 1000)
##    cost misclass
## 1 1e-04      480
## 2 1e-03       86
## 3 1e-02       92
## 4 1e+00       89
## 5 5e+00       92
## 6 1e+01       92
## 7 1e+02       89
## 8 1e+03       89
## 9 1e+04       91

Podemos ver que con un costo de 0.01 se clasificaron mal 85 observaciones.

set.seed(5)
x1_test=runif (1000) -0.5 
x2_test=runif (1000) -0.5

temp_test <- x1_test-x2_test

y_test <- rep(0, 1000)

for (i in 1:1000) {
  if (temp_test[i] > 0.1) {
    y_test[i] <- 1
  }
  else if (temp_test[i] < -0.1) {
    y_test[i] <- 0
  }
  else {
    y_test[i] <- runif(1) > 0.5
  }
}
  
df3 <- data.frame(x1_test,x2_test,y_test)

df3$y_test <- as.factor(df3$y_test)

ggplot(df3, aes(x=x1_test,y=x2_test,col=y_test)) +
  geom_point()

names(df3) <- c("x1","x2","y")
costs <- c(0.0001, 0.001, 0.01, 1, 5, 10, 100, 1000, 10000)
test.err <- rep(NA, length(costs))
for (i in 1:length(costs)) {
    svm.fit <- svm(y ~ ., data = df2, kernel = "linear", cost = costs[i])
    pred <- predict(svm.fit, df3)
    test.err[i] <- sum(pred != df3$y)
}
data.frame(cost = costs, misclass = test.err)
##    cost misclass
## 1 1e-04      486
## 2 1e-03       91
## 3 1e-02       92
## 4 1e+00       95
## 5 5e+00       94
## 6 1e+01       96
## 7 1e+02       96
## 8 1e+03       96
## 9 1e+04       96

Podemos ver que en este caso el mejor costo fue dado por \(c=1\), con 92 errores a la hora de predicción

Podemos apreciar que aunque el costo en ambos casos fue diferente, realmente no difiere demasiado, ya que a partir de un costo de 0.001 no se presentan grandes cambios en los errores,

#### 7

Auto <- Auto[,-10] %>%
    mutate(highmpg = as.integer(mpg > median(mpg))) %>%
    mutate(highmpg = factor(highmpg))
head(Auto)
##   mpg cylinders displacement horsepower weight acceleration year origin
## 1  18         8          307        130   3504         12.0   70      1
## 2  15         8          350        165   3693         11.5   70      1
## 3  18         8          318        150   3436         11.0   70      1
## 4  16         8          304        150   3433         12.0   70      1
## 5  17         8          302        140   3449         10.5   70      1
## 6  15         8          429        198   4341         10.0   70      1
##                        name highmpg
## 1 chevrolet chevelle malibu       0
## 2         buick skylark 320       0
## 3        plymouth satellite       0
## 4             amc rebel sst       0
## 5               ford torino       0
## 6          ford galaxie 500       0
svml <- tune(svm, highmpg ~ ., data = Auto, kernel = "linear", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000)))
summary(svml)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     1
## 
## - best performance: 0.01275641 
## 
## - Detailed performance results:
##    cost      error dispersion
## 1 1e-03 0.09205128 0.05715898
## 2 1e-02 0.07416667 0.04448352
## 3 1e-01 0.04865385 0.03720214
## 4 1e+00 0.01275641 0.01808165
## 5 5e+00 0.02044872 0.02020886
## 6 1e+01 0.02551282 0.01709615
## 7 1e+02 0.03826923 0.01345323
## 8 1e+03 0.03826923 0.01345323

Se puede apreciar como con un costo de 1 se tiene el menor error.

  • Kernel polinomial
svmk <- tune(svm, highmpg ~ ., data = Auto, kernel = "polynomial", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000), degree = c(2, 3, 4)))
summary(svmk)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost degree
##  1000      2
## 
## - best performance: 0.2472436 
## 
## - Detailed performance results:
##     cost degree     error dispersion
## 1  1e-03      2 0.5635256 0.03971883
## 2  1e-02      2 0.5635256 0.03971883
## 3  1e-01      2 0.5635256 0.03971883
## 4  1e+00      2 0.5635256 0.03971883
## 5  5e+00      2 0.5635256 0.03971883
## 6  1e+01      2 0.5251282 0.11987686
## 7  1e+02      2 0.2959615 0.06069283
## 8  1e+03      2 0.2472436 0.05930020
## 9  1e-03      3 0.5635256 0.03971883
## 10 1e-02      3 0.5635256 0.03971883
## 11 1e-01      3 0.5635256 0.03971883
## 12 1e+00      3 0.5635256 0.03971883
## 13 5e+00      3 0.5635256 0.03971883
## 14 1e+01      3 0.5635256 0.03971883
## 15 1e+02      3 0.3341667 0.08484496
## 16 1e+03      3 0.2551282 0.06145828
## 17 1e-03      4 0.5635256 0.03971883
## 18 1e-02      4 0.5635256 0.03971883
## 19 1e-01      4 0.5635256 0.03971883
## 20 1e+00      4 0.5635256 0.03971883
## 21 5e+00      4 0.5635256 0.03971883
## 22 1e+01      4 0.5635256 0.03971883
## 23 1e+02      4 0.5635256 0.03971883
## 24 1e+03      4 0.5302564 0.09378263

Aquí vemos que el menor error se da con un costo de 1000 de grado 2

  • Kernel radial
svmr <- tune(svm, highmpg ~ ., data = Auto, kernel = "radial", ranges = list(cost = c(0.001,0.01, 0.1, 1, 5, 10, 100, 1000), gamma = c(0.01, 0.1, 1, 5, 10, 100)))
summary(svmr)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost gamma
##   100  0.01
## 
## - best performance: 0.01532051 
## 
## - Detailed performance results:
##     cost gamma      error dispersion
## 1  1e-03 1e-02 0.59147436 0.06357909
## 2  1e-02 1e-02 0.59147436 0.06357909
## 3  1e-01 1e-02 0.08935897 0.03875505
## 4  1e+00 1e-02 0.07647436 0.03804385
## 5  5e+00 1e-02 0.05608974 0.03963690
## 6  1e+01 1e-02 0.03314103 0.02424635
## 7  1e+02 1e-02 0.01532051 0.01788871
## 8  1e+03 1e-02 0.03051282 0.02878864
## 9  1e-03 1e-01 0.59147436 0.06357909
## 10 1e-02 1e-01 0.24179487 0.10616790
## 11 1e-01 1e-01 0.08166667 0.03582639
## 12 1e+00 1e-01 0.05608974 0.03963690
## 13 5e+00 1e-01 0.02801282 0.02509397
## 14 1e+01 1e-01 0.02288462 0.02505026
## 15 1e+02 1e-01 0.03307692 0.03404377
## 16 1e+03 1e-01 0.03307692 0.03404377
## 17 1e-03 1e+00 0.59147436 0.06357909
## 18 1e-02 1e+00 0.59147436 0.06357909
## 19 1e-01 1e+00 0.59147436 0.06357909
## 20 1e+00 1e+00 0.06365385 0.04003438
## 21 5e+00 1e+00 0.06615385 0.04478910
## 22 1e+01 1e+00 0.06615385 0.04478910
## 23 1e+02 1e+00 0.06615385 0.04478910
## 24 1e+03 1e+00 0.06615385 0.04478910
## 25 1e-03 5e+00 0.59147436 0.06357909
## 26 1e-02 5e+00 0.59147436 0.06357909
## 27 1e-01 5e+00 0.59147436 0.06357909
## 28 1e+00 5e+00 0.52775641 0.07143278
## 29 5e+00 5e+00 0.52262821 0.07036637
## 30 1e+01 5e+00 0.52262821 0.07036637
## 31 1e+02 5e+00 0.52262821 0.07036637
## 32 1e+03 5e+00 0.52262821 0.07036637
## 33 1e-03 1e+01 0.59147436 0.06357909
## 34 1e-02 1e+01 0.59147436 0.06357909
## 35 1e-01 1e+01 0.59147436 0.06357909
## 36 1e+00 1e+01 0.55583333 0.05641062
## 37 5e+00 1e+01 0.54820513 0.06563816
## 38 1e+01 1e+01 0.54820513 0.06563816
## 39 1e+02 1e+01 0.54820513 0.06563816
## 40 1e+03 1e+01 0.54820513 0.06563816
## 41 1e-03 1e+02 0.59147436 0.06357909
## 42 1e-02 1e+02 0.59147436 0.06357909
## 43 1e-01 1e+02 0.59147436 0.06357909
## 44 1e+00 1e+02 0.59147436 0.06357909
## 45 5e+00 1e+02 0.59147436 0.06357909
## 46 1e+01 1e+02 0.59147436 0.06357909
## 47 1e+02 1e+02 0.59147436 0.06357909
## 48 1e+03 1e+02 0.59147436 0.06357909

Con un costo de 100 y un gama de 0.01 se obtuvo el menor error en el modelo

svm.linear <- svm(highmpg ~ ., data = Auto, kernel = "linear", cost = 1)
svm.poly <- svm(highmpg ~ ., data = Auto, kernel = "polynomial", cost = 100, degree = 2)
svm.radial <- svm(highmpg ~ ., data = Auto, kernel = "radial", cost = 100, gamma = 0.01)
  • Lineal

  • Polinomial

  • Radial

8

data('OJ')

inTrain <- sample(nrow(OJ), 800, replace = FALSE)

training <- OJ[inTrain,]
testing <- OJ[-inTrain,]
svm_linear <- svm(Purchase ~ ., data = training,
                  kernel = 'linear',
                  cost = 0.01)
summary(svm_linear)
## 
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "linear", 
##     cost = 0.01)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  0.01 
## 
## Number of Support Vectors:  445
## 
##  ( 222 223 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM

Podemos apreciar que se tienen 440 vectores de soporte de los cuales 219 pertenecen a CH y 221 a MM

train.pred <- predict(svm_linear, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 411  76
##         MM  64 249
## [1] "Accuracy: 0.825 Sensitivity: 0.8653 Specificity: 0.7662"
# Calculando el error de entrenamiento
a <- 1-matriz$overall[1]
names(a) <- "error"
a
## error 
## 0.175
test.pred <- predict(svm_linear, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 158  27
##         MM  20  65
## [1] "Accuracy: 0.8259 Sensitivity: 0.8876 Specificity: 0.7065"
# Calculando el error del set de prueba
a <- 1-matriz$overall[1]
names(a) <- "error"
a
##     error 
## 0.1740741

Podemos ver que en ambos casos el error es del 17%

to = tune(svm, Purchase ~ ., data = training, kernel = "linear", ranges = list(cost = 10^seq(-2, 
    1, by = 0.25)))
summary(to)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##      cost
##  1.778279
## 
## - best performance: 0.16625 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.18000 0.04495368
## 2   0.01778279 0.17500 0.03864008
## 3   0.03162278 0.17250 0.04518481
## 4   0.05623413 0.17000 0.04297932
## 5   0.10000000 0.17250 0.04440971
## 6   0.17782794 0.17125 0.04566256
## 7   0.31622777 0.17625 0.04619178
## 8   0.56234133 0.17125 0.04450733
## 9   1.00000000 0.17000 0.04721405
## 10  1.77827941 0.16625 0.04678927
## 11  3.16227766 0.16750 0.04937104
## 12  5.62341325 0.16750 0.05075814
## 13 10.00000000 0.16875 0.05008673
svm.linear <- svm(Purchase ~ ., kernel = "linear", data = training, cost = to$best.parameters$cost)
train.pred = predict(svm.linear, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 415  72
##         MM  60 253
## [1] "Accuracy: 0.835 Sensitivity: 0.8737 Specificity: 0.7785"
## error 
## 0.165
test.pred <- predict(svm_linear, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 158  27
##         MM  20  65
## [1] "Accuracy: 0.8259 Sensitivity: 0.8876 Specificity: 0.7065"
##     error 
## 0.1740741

Podemos ver que el error de entrenamiento baja un poco, pero el de prueba se mantiene

svm.radial = svm(Purchase ~ ., data = training, kernel = "radial")
summary(svm.radial)
## 
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  372
## 
##  ( 190 182 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM
train.pred = predict(svm.radial, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 428  78
##         MM  47 247
## [1] "Accuracy: 0.8438 Sensitivity: 0.9011 Specificity: 0.76"
##   error 
## 0.15625
test.pred = predict(svm.radial, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 159  34
##         MM  19  58
## [1] "Accuracy: 0.8037 Sensitivity: 0.8933 Specificity: 0.6304"
##     error 
## 0.1962963
to = tune(svm, Purchase ~ ., data = training, kernel = "radial", ranges = list(cost = 10^seq(-2, 
    1, by = 0.25)))
summary(to)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##      cost
##  3.162278
## 
## - best performance: 0.17625 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.40625 0.04903584
## 2   0.01778279 0.40625 0.04903584
## 3   0.03162278 0.32625 0.06440163
## 4   0.05623413 0.20250 0.05263871
## 5   0.10000000 0.18500 0.04281744
## 6   0.17782794 0.18375 0.04291869
## 7   0.31622777 0.17750 0.04322101
## 8   0.56234133 0.17750 0.05027701
## 9   1.00000000 0.17875 0.05369991
## 10  1.77827941 0.17750 0.05163978
## 11  3.16227766 0.17625 0.05350558
## 12  5.62341325 0.17625 0.05015601
## 13 10.00000000 0.17750 0.05130248
svm.radial = svm(Purchase ~ ., data = training, kernel = "radial", cost = to$best.parameters$cost)
train.pred = predict(svm.radial, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 430  69
##         MM  45 256
## [1] "Accuracy: 0.8575 Sensitivity: 0.9053 Specificity: 0.7877"
##  error 
## 0.1425
test.pred = predict(svm.radial, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 160  35
##         MM  18  57
## [1] "Accuracy: 0.8037 Sensitivity: 0.8989 Specificity: 0.6196"
##     error 
## 0.1962963
svm.poly = svm(Purchase ~ ., data = training, kernel = "poly", degree = 2)
summary(svm.poly)
## 
## Call:
## svm(formula = Purchase ~ ., data = training, kernel = "poly", 
##     degree = 2)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  polynomial 
##        cost:  1 
##      degree:  2 
##      coef.0:  0 
## 
## Number of Support Vectors:  459
## 
##  ( 235 224 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM
train.pred = predict(svm.poly, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 430 107
##         MM  45 218
## [1] "Accuracy: 0.81 Sensitivity: 0.9053 Specificity: 0.6708"
## error 
##  0.19
test.pred = predict(svm.poly, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 171  46
##         MM   7  46
## [1] "Accuracy: 0.8037 Sensitivity: 0.9607 Specificity: 0.5"
##     error 
## 0.1962963
to = tune(svm, Purchase ~ ., data = training, kernel = "poly", degree = 2, 
    ranges = list(cost = 10^seq(-2, 1, by = 0.25)))
summary(to)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##    10
## 
## - best performance: 0.185 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.40375 0.05894029
## 2   0.01778279 0.37750 0.06449591
## 3   0.03162278 0.36875 0.06594663
## 4   0.05623413 0.34625 0.06237487
## 5   0.10000000 0.31750 0.06487167
## 6   0.17782794 0.24125 0.06292908
## 7   0.31622777 0.21750 0.04609772
## 8   0.56234133 0.21125 0.04619178
## 9   1.00000000 0.20875 0.04641674
## 10  1.77827941 0.19375 0.04649149
## 11  3.16227766 0.19000 0.03809710
## 12  5.62341325 0.18750 0.04124790
## 13 10.00000000 0.18500 0.03855011
svm.poly = svm(Purchase ~ ., data = training, kernel = "poly", degree = 2, cost = to$best.parameters$cost)
train.pred = predict(svm.poly, training)
matriz <- confusionMatrix(train.pred, training$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 433  74
##         MM  42 251
## [1] "Accuracy: 0.855 Sensitivity: 0.9116 Specificity: 0.7723"
## error 
## 0.145
test.pred = predict(svm.poly, testing)
matriz <- confusionMatrix(test.pred, testing$Purchase)
matriz$table
##           Reference
## Prediction  CH  MM
##         CH 164  37
##         MM  14  55
## [1] "Accuracy: 0.8111 Sensitivity: 0.9213 Specificity: 0.5978"
##     error 
## 0.1888889

Podemos ver como el polinomial usando el mejor costo es el que ofrece un mejor resultado de entrenamiento pero el lineal es el que tiene menor error de prueba

Ensayo


Medellín es una ciudad la cual debido a su crecimiento poblacional y posición geográfica, presenta en general unos altos índices de contaminación, en particular, en algunas épocas del año esta crisis ambiental se agudiza. Más del 80% de las personas que viven en áreas urbanas que monitorean la contaminación del aire están expuestas a niveles de calidad del aire que exceden los límites de la Organización Mundial de la Salud (OMS, 2016). Esto hace necesario que se busquen mecanismos que ayuden a mejorar la problemática de la calidad del aire.

Por un lado tenemos que debido a la topografía de Medellín, y en general del valle de Aburrá, en épocas del año en donde las condiciones climáticas son propicias para que el aire no pueda calentarse lo suficiente, y en consecuencia no pueda ascender para que esta sea arrastrado fuera del Valle por los vientos alisios. Desde este punto de vista, es importante profundizar en la investigación de los modelos predictivos meteorológicos, los cuales, al poder darnos alertas mucho más precisas, podrían usarse para tomar medidas oportunas y no especulativas de cara a afrontar el problema.

Para llevar esto a cabo se puede tomar como referencia al trabajo realizado por Y. Zheng et al. en el año 2015, en donde pronosticaron la lectura de una estación de monitoreo de la calidad del aire de las siguientes 48 horas utilizando un modelo basado en datos (data-driven). Este sistema cuenta con un predictor de regresión lineal para modelar los factores locales de la calidad del aire. Para modelar los factores globales se usa una red neuronal, el cual es un predictor espacial. También cuenta con un agregador dinámico el cual combina las predicciones de los predictores teniendo en cuenta también los datos metereológicos y por último cuenta con un predictor para capturar los cambios repentinos en la calidad del aire.

Pero para esto tal y como subrayan en el año 2017 Kaur et al. existe una necesidad por mejorar la toma de información, ya que en definitiva, de esta dependen todos los modelos que se planteen. En este apartado es importante que Medellín cuente con un sistema de medición en tiempo real en donde no solo se monitoreen algunas áreas de la ciudad, si no donde podamos tener un entendimiento de que es lo que está ocurriendo realmente en las diferentes calles de Medellín, para esto podríamos usar diferentes dispositivos ayudándonos del internet de las cosas. Uno de estos ejemplos lo presenta Wooden en el año 2019, en donde nos muestra como Intel y Bosch desarrollaron un sistema al que denominaron The Intel®-based Bosch Air Quality Micro Climate Monitoring System (MCMS), el cual sin ocupar mucho espacio, podría usarse en diferentes zonas de la ciudad, para saber no solamente, si en general en Medellín hay mucha contaminación, si no en qué lugares específicamente es donde más contaminación se está presentando, en que calles en específico y a que horas del día, para poder tomar más medidas en estos lugares, ya que la contaminación que allí se genera termina teniendo un impacto en una zona mucho más amplia y afecta directamente a una gran parte de la población. Esto no solo ayudaría a tomar medidas, si no que permitiría evaluarlas para saber si realmente estas eran el problema principal.

  • Referencias:

Organización Mundial de la Salud. (2016). Ambient Air Pollution: A global assessment of exposure and burden of disease. World Health Organization, 1–131. Retrieved from www.who.int.org

Y. Zheng, X. Yi, M. Li, R. Li, Z. Shan, E. Chang and T. Li, “Forecasting Fine-Grained Air Quality Based on Big Data,” Proceeding KDD ’15 Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 2267-2276, August 10, 2015.

Kaur, Gaganjot & Gao, Jerry & Chiao, Sen & Lu, Shengqiang & Xie, Gang. (2017). Air Quality Prediction: Big data and Machine Learning Approaches.

Wooden, A. (2019). Tackling air pollution with big data analytics. Intel. Recuperado de https://www.intel.co.uk